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Abstract 

We present dynamic Monte Carlo simulations of a lattice-gas model for bromine electrode- 
position on single-crystal silver (100). This system undergoes a continuous phase transition be- 
tween a disordered phase at low electrode potentials and a commensurate c(2 x 2) phase at high 
potentials. The lattice-gas parameters are determined by fitting simulated equilibrium adsorp- 
tion isotherms to chronocoulometric data, and free-energy barriers for adsorption/desorption 
and lateral diffusion are estimated from ab initio data in the literature. Cyclic voltammo- 
grams in the quasi-static limit are obtained by equilibrium Monte Carlo simulations, while for 
nonzero potential scan rates we use dynamic Monte Carlo simulation. The butterfly shapes of 
the simulated voltammograms are in good agreement with experiments. Simulated potential- 
step experiments give results for the time evolution of the Br coverage, as well as the c(2 x 2) 
order parameter and its correlation length. During phase ordering following a positive potential 
step, the system obeys dynamic scaling. The disordering following a negative potential step is 
well described by random desorption with diffusion. Both ordering and disordering processes 
are strongly influenced by the ratio of the time scales for desorption and diffusion. Our re- 
sults should be testable by experiments, in particular cyclic voltammetry and surface X-ray 
scattering. 

Keywords: Bromine adsorption; Continuous phase transition; Cyclic Voltammetry; Dynamic 
Monte Carlo simulation; Lattice-gas model; Potential-step experiments; 

1 Introduction 

The electrodeposition of Br on single-crystal Ag(100) from aqueous solution is a simple example 
of anion adsorption which has been extensively studied, both by classical electrochemical methods 
jj], ||, ||, ||, H and by techniques such as in situ surface X-ray scattering (SXS) jf, H §| an d X-ray 
absorption fine structure (XAFS) |Q. 

Cyclic voltammetry (CV) shows a typical butterfly structure with a broad pre-wave in the 
negative-potential region and a sharp peak at more positive potentials Jj], ||, |[ 0], The recent 
SXS experiments by Ocko, Wang, and Wandlowski [|| showed that the sharp peak corresponds to 
a continuous phase transition in the layer of adsorbed Br. At this transition the adlayer changes 
its structure from a disordered two-dimensional "gas" on the negative-potential side, to an ordered, 
commensurate c(2 x 2) phase with a Br coverage of 1/2 monolayer on the positive side. The pre-wave 
lies wholly in the disordered-phase region, and it was previously suggested that it was caused by in- 
teractions with surface water ffl. However, equilibrium Monte Carlo (MC) simulations of water-free 
lattice-gas models of the Br adlayer by Koper || [7), as well as our own work reported here, produce 
quasi-equilibrium cyclic voltammograms (CVs) of essentially the same shape as the experiments, 
indicating that surface water is not needed to reproduce the pre-wave. Rather, we believe that 
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the pre- wave is due to short-range correlations in the disordered adlayer [g, ||. In this system, the 
structures that give rise to these correlations locally resemble low-temperature ordered phases with 
coverage 1/4. 

The simplicity of the equilibrium properties of this system makes it a prime candidate for dynamic 
studies. Here we present results from such a study by dynamic MC simulation |)| of a lattice-gas 
model. This method has the advantage over mean-field rate-equation approaches that it properly 
accounts for the effects of local fluctuations in the adlayer structure. Further details and additional 
results will be presented elsewhere . Animations and additional figures are available on the World 
Wide Web Q. 

The rest of this paper is organized as follows. The lattice-gas model is presented in Sec. ||, followed 
by equilibrium simulations which are used to estimate the lattice-gas model parameters by fitting to 
adsorption isotherms from chronocoulometry experiments || ||[ ||. In Sec. || we present simulated 
CVs (Sec. 3J) and potential-step experiments (Sec. O). From the the potential-step simulations 
we provide current transients, which are easy to measure in electrochemical experiments, as well as 
predictions for the time evolution of the intensity of the SXS scattering peak corresponding to the 
c(2 x 2) phase. Although dynamic SXS data have yet to be obtained for this system, they were 
obtained for others iQ, and we believe it is only a matter of time before such experimental results 
become available. 



2 Model and Equilibrium Results 

We use a lattice-gas model similar to the one used by Koper j6). It consists of an L x L square 
array of Br adsorption sites, corresponding to the four-fold hollow sites on the Ag(100) surface. The 
configurational energy of the Br adlayer is given by the grand-canonical lattice-gas Hamiltonian, 

L 2 

H = -y j fajCjCj -jLjTci . (l) 

Here i and j denote adsorption sites, Cj is the occupation at site i, which is either (empty) or 1 
(occupied), J2i<j ls a sum over an pairs of sites, <f>ij is the lateral interaction energy of the pair 
(«, j), and ~p is the electrochemical potential. The sign conventions are such that 4>ij < denotes a 
repulsive interaction, and Jl > favors adsorption. We measure fej and Jl in units of meV/pair and 
meV/particle, respectively. (For brevity, both will be written simply as meV.) To reduce finite-size 
effects, we use periodic boundary conditions. 

In the weak-solution approximation, Jl is related to the electrode potential by 

[CI 

JJ = JJ + fc B Tln -—r - ejE , (2) 

where ~p is an arbitrary reference level, fee is Boltzmann's constant, T is the absolute temperature, 
e is the elementary charge unit, [C] is the concentration of Br~ in solution, [Co] is an arbitrary 
reference concentration, 7 is the electrosorption valency, and E is the electrode potential in mV. 

Previously Koper has explored the effects of finite nearest-neighbor repulsion and screened dipole- 
dipole interactions in this model J^]. As his results indicate only minor effects of the finite nearest- 
neighbor repulsion and screening, we here use a simplified model with nearest-neighbor exclusion 
and unscreened dipole-dipole interactions. Thus, the long-range part of the interaction energy is 
<j>{r) = 2 3 / 2 </> nnn r~ 3 for r > \/2, where r is the separation of an interacting Br pair, measured in 
units of the Ag(100) lattice spacing (a = 2.889 A [§), and 4>nnn (which is negative) is the lateral 
dipole-dipole repulsion between next-nearest neighbors. For computational convenience we cut off 
the long-range interaction for r > 5. This underestimates the total interaction energy of a fully 
occupied c(2 x 2) layer by about 13%. 

The nearest-neighbor exclusions give rise to two different sublattices of possible adsorption sites 
(like the black and white squares on a chessboard). Each sublattice (labeled A and B, respectively) 
corresponds to one of two degenerate c(2 x 2) phases. The A sublattice coverage is the fraction of 
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occupied sites on sublattice A, denned as 6a = A^ 1 YlieA Ci ' wnere is the number of sites on 
sublattice A and YlieA runs over a H s it es on the sublattice. We define 6 b analogously. 

The sublattice coverages combine to give two observables of interest: the total Br coverage 

6 = (6a + Bb)/2, and the "staggered" coverage 8s = 6a — 6b, which is the order parameter for 
the c(2 x 2) phase. While 6 can be experimentally obtained by standard electrochemical methods, 
as well as from the integer-order peaks in scattering data, 6s is proportional to the square root of 
the intensity in the half-order diffraction peaks that correspond to the c(2 x 2) phase [Q. 

To estimate the parameters in the lattice-gas model, <p nnn and 7, we performed standard equi- 
librium MC simulations jl3|] at room temperature (k^T = 25 meV or T w 290 K) to obtain 6(7*) 
for different parameter values. These simulated isotherms were then compared with experimental 
chronocoulometry data for three different electrolyte concentrations M, [jj, and the best parameter 
values were determined by a nonlinear fit The resulting values are </> nnn = — 26 ± 2 meV and 

7 = — 0.73± 0.03. These results are consistent with those found previously [I], ||. The fitted MC 
isotherm is shown together with the experimental isotherms in Fig. [jj. Especially for the two lowest 
concentrations, the agreement is excellent over the whole range of electrode potentials. The kinks in 
the isotherms, observed at 6 C « 0.37 for both the experiments and simulations, correspond to the 
order-disorder phase transition. To within our statistical uncertainty, this value of 6 C is consistent 
with previous results for a s qua re-lattice model with only nearest-neighbor exclusion, known as the 
"hard-square model" Q |f @. 

The simulation data for 6(/2) and &s(p) with the best-fit interaction parameter (f> nim = —26 meV 
are shown together in Fig. |[ As p approaches its critical value at p c — 180 ±5 meV from the positive 
side, 6s oc (p — pc) 1 / 8 - This behavior of the order parameter is consistent with the Ising universality 
class, to which this system belongs, and was confirmed by SXS experiments ||. Typical equilibrium 
configurations are shown as insets in Fig. |^ for the disordered phase at p = +100 meV and for the 
ordered phase at p — +400 meV. 

A simulated quasi-equilibrium CV, corresponding to a vanishing potential-sweep rate, can be 
obtained from the MC equilibrium isotherm as 

J - jeM dQ - 7 eM de ^ dE - 1 2 e 2 M dQ dE oc d ® (3) 
dt dp dE dt dp dt dp ' 

where J is the voltammetric current density (oxidation currents positive), M is the total number 
of adsorption sites per unit area {aT 2 = 1.198 x 10 15 sites/cm 2 ), and dE/dt is the sweep rate. 
This limiting CV is shown in Figure |^(a) together with CVs for nonzero sweep rates, which are 



obtained by dynamic MC simulations discussed in Sec. 3.1. It exhibits the same broad pre-wave and 
sharp peak seen in experiments [|, |j| [J, The broad prewave in the simulated CV is caused by 
configurational fluctuations in the disordered phase, which locally resemble low-temperature ordered 
phases [jlT], |l8|], in particular p(2 x 2) and c(4 x 2) with 6 = 1/4 0. Such fluctuations of local 
short-range order in the disordered phase are clearly seen in the inset in Fig. || for p = +100meV, and 
they lead to visible anisotropy in simulated diffuse SXS scattering intensities QlCJ] . This anisotropy 
should be experimentally observable as well. 



3 Dynamic Simulations 

The dynamics of the Br adsorption and desorption processes under C V and potential-step conditions 
were simulated with a dynamic MC algorithm including adsorption and desorption events, as well 
as nearest- and next-nearest-neighbor lateral diffusion of the adsorbed Br. Bulk diffusion in the 
solution is neglected, corresponding to a well-stirred system. Each such single microscopic move, 
which we label by the index A, connects an initial lattice-gas state, /, to a final state, F. The energies 
of these states, Ui and Uf, are obtained by applying the lattice-gas Hamiltonian, Eq. ([[]), to the 
corresponding configurations. An intermediate transition state of higher energy, T\, is associated 
with the move A. This intermediate state cannot be represented by a lattice-gas configuration, and we 
associate with it a "bare" free-energy barrier, A A . Using a symmetric Butler- Volmer approximation 
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H |l(| |l9|], we can then approximate the free energy of the transition state as 



^ = ^^+A A . (4) 

Although other choices of the transition probability are also found in the literature |m| , we here 
approximate the probability R(F\I) of making a transition from I to F during a single MC time 
step by the one-step Arrhenius rate || ||l], ||^] 

R{m . „_ . „ ffip » p (_^) , 

where the dimensionless rate constant v relates the simulation time scale, measured in MC steps per 
site (MCSS), to the experimental time scale in seconds. In our simulations we used v = 1. 

We performed the simulation with a simple discrete-time dynamical algorithm, in which new 
configurations are randomly chosen from a weighted list of microscopic moves. First, a particular 
lattice site is chosen at random. The configuration can then be changed only through moves which 
include the chosen site. If the site is empty, adsorption is attempted and accepted with probability 
given by Eq. (||). If any of the nearest-neighbor sites are occupied, this acceptance probability is 
zero. If the chosen site is occupied, only desorption or lateral diffusion may be attempted. A list is 
kept of all the possible moves, and one of them is chosen according to the corresponding acceptance 
probabilities. The probability to remain in the initial state is R(I\I) = 1 — J^f=£t Further 
details on the simulation algorithm and its implementation are given elsewhere [g[ |l0[ . 

While this algorithm is of limited accuracy for very fast processes and requires an unnecessarily 
large amount of computer time for very slow processes J2l], |23|, |J, ^5[ |(| |27|] , it has the advantages 
that it is easy to program and is readily adapted to simulations in which parameters (such as p,) 
change with time ^7|. It is thus well suited for dynamic CV simulations. 

The free-energy barriers associated with the different simulation moves are A nn for nearest- 
neighbor diffusion, A nnn for next-nearest neighbor diffusion, and A a for adsorption/desorption. The 
rough estimates used here, A nn = 100 meV and A nnn = 200 meV, are based on ab-initio calculations 
of binding energies for a single Br ion on a Ag(100) substrate in vacuum |28). The difference in 
binding energy between the bridge site and the four-fold hollow site gives A nn , and the difference 
between the on-top site and the four-fold hollow site gives A nnn . Theoretical estimates of A a 
are extremely sensitive to the ion-surface and ion- water interactions pg| |. Calculated potentials of 
mean force for halide ions in water near a Cu(100) surface [^9| indicate that values between 200 
and 500 meV are not unreasonable. To optimize the simulation speed, we chose A a = 300 meV. 
This value is as low as possible while remaining consistent with our expectation that it should be 
significantly larger than A nnn . 



3.1 Simulated Cyclic Voltammograms 

CV experiments were simulated on systems with L = 256 and 128 by first equilibrating at p, = 
—400 meV and then ramping ft linearly in time up to +600 meV, and back down to —400 meV. 
Simulated CV currents for sweep rates between 3xl0 -4 and 3xl0 -3 meV/MCSS, divided by the 
sweep rate in order to be easily compared on the same scale, are shown in Fig. ||(a). The limiting 
curve for vanishing sweep rate is given by the quasi-equilibrium CV, obtained from the equilibrium 
simulations described in Sec. ||. A considerable asymmetry in the peak shape for the positive-going 
and negative-going scans is notable. 

While the pre-waves are relatively little affected by the sweep rate, the peak corresponding to 
the phase transition is significantly lowered and shifted in the scan direction. The dependence of 
the peak separation on the sweep rate, which is a form of hysteresis J30), is shown in Fig. ||(b). 

By comparing the scan-rate dependence of the separation between the positive-going and negative- 
going peak positions in simulations and experiments, one can in principle obtain a rough estimate of 
the relation between simulation time and physical time, provided the free-energy barriers used in the 
simulations are in a reasonably realistic proportion to each other. The peak separations observed in 
our simulations, Fig. ^(b), are large compared to typical experimental values of a few tens of meV 



4 



for scan rates in the range 1-10 mV/s. This most likely indicates that even our slowest simulated 
scan rates correspond to faster scans in the real system. Direct comparisons between simulation and 
experimental dynamic effects are left for future research. 



3.2 Simulated Potential Steps 

Dynamic potential-step simulations were performed on systems with L = 256. They began by 
equilibrating the system at room temperature and potential After equilibration, £l was instanta- 
neously stepped to a new value, J12, and the simulation was continued by the dynamic MC algorithm 
described above. To illustrate the dynamics of the phase ordering and disordering processes, we here 
show results for two different potential steps: one from the disordered phase into the ordered phase, 
and one from the ordered phase into the disordered phase. Additional potential-step simulations 
and corresponding time-dependent SXS scattering intensities will be reported elsewhere [ |To| , 



3.2.1 Disorder-to-order step 

For the disorder-to-order step we used jli = —200 meV and /I2 = +600 meV. At p,i the coverage is 
close to zero, while p,2 is far into the ordered phase, approximately 420 meV past the phase transition 
at /2 C . In Fig. || we show both 8 and 9s vs time. With deep steps like this one, the desorption 
rate is negligible. Thus, the adsorption dynamics are essentially described by the random sequential 
adsorption with diffusion (RSAD) process jH], |32). The coverage quickly reaches the jamming 



coverage for random sequential adsorption of hard squares without diffusion, 8j = 0.364 |33| |3j 
(only slightly below the critical coverage, 9 C ~ 0.37). At this coverage both sublattices contain 
small, uncorrelated domains of the two degenerate ordered phases, separated by domain walls. The 
domain walls consist of empty sites, most of which cannot be filled, due to the nearest-neighbor 
exclusion. 

At later times, the coverage can increase only where domain walls move together, opening a gap 
large enough to fit an additional Br. The domain-wall motion responsible for opening additional 
adsorption sites proceeds almost exclusively through nearest-neighbor diffusion between adjacent 
domains. The coverage approaches the equilibrium value, 9 cq = 1/2, and the total length of 
interfaces decreases. The order-parameter correlation length, D, is experimentally measurable as 
the inverse width of the half-order diffuse SXS scattering maxima |g, [ijj . Since D is also proportional 
to the inverse of the interfacial length per unit area |3^], it can in the present case be estimated 
from the coverage as D = (9 cq — 9) _1 |n], ^|. The dynamic scaling theory of the kinetics of phase 
transitions predicts that, for systems with nonconserved order parameter (such as the one considered 
here) undergoing phase ordering, D should grow with time as t 1 / 2 |56|. The inset in Fig. ^ shows 
D as a function of t 1 ^ 2 . The agreement with the expected dynamic scaling behavior is excellent, as 
has also been found in previous simulations of RSAD J3l], |32^| . 

For an infinite system, D would continue to grow without bound; however, when D reaches the 
same order of magnitude as L, the system enters the phase-selection regime in which one of the 
two degenerate ordered phases grows to fill the whole system. Our step simulations were too short 
to study the phase-selection regime in detail. Results for shallower positive-going potential steps, 
which are somewhat less clear-cut than for the deep step shown here, will be reported elsewhere [[10||. 



3.2.2 Order-to-disorder step 

For deep order-to-disorder steps, the disordering process is rather simple. Particles desorb at a 
roughly constant rate, leading to an essentially exponential relaxation to equilibrium |jlo|| . 

After a shallow order-to-disorder step, fli = +600 meV and ^2 = +100 meV, the behavior is 
more interesting. The system starts with all particles on sublattice A and relaxes to a disordered 
phase with 9 w 1/4, as shown in Fig. |[ There are four different dynamical regimes. In the first 
regime, particles simply desorb from sublattice A so that d&s/dt f=s 2dO/dt. As more sites become 
vacant, small domains are formed on sublattice B both by lateral diffusion from sublattice A and 
by adsorption from the solution. In this second regime, both desorption and diffusion contribute 
significantly to the disordering process. In the third regime, the adsorption and desorption rates are 
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almost equal and relatively slow, and diffusion is the dominant contribution to the disordering. Here 
dQs/dt becomes significantly larger than d<d/dt, until the two sublattices become approximately 
equally populated in the fourth regime, yielding 8s ~ 0. 



4 Conclusion 

In this brief paper we have presented results from dynamic MC simulations of CV and potential-step 
experiments for Br electrosorption on Ag(100) single-crystal electrodes. The simulations show the 
complicated interplay between adsorption, desorption, and lateral diffusion even in this relatively 
simple electrochemical system. 

The dynamic MC algorithm requires estimates of free-energy barriers for adsorption/desorption 
and lateral diffusion. While such barrier estimates constitute the most uncertain part of a dy- 
namic MC simulation, we were able to obtain reasonable numbers from ab initio calculations in the 
literature [|§ |§. 

As the order parameter 0s and its correlation length are measurable in SXS experiments 
||, ||, H^|, the dynamical phenomena predicted by our simulations should be observable in future 
experiments. Further details and simulations of time-dependent diffuse SXS scattering intensities 
will be reported elsewhere [{Tcfl . 
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Figure 1: A single equilibrium Monte Carlo (MC) simulation isotherm (L = 32, T ss 290 K), 
fit simultaneously to coverage isotherms at three different Br~ concentrations, obtained by chrono- 
coulometry |3[ ||. The parameters, 4> nxm — —26 ± 2 meV and 7 = —0.73 ± 0.03, were obtained by a 
nonlinear fit to the experimental data. 
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Figure 2: Equilibrium Monte Carlo isotherms for L = 32, T w 290 K, and n nn = —26 meV. 
A continuous phase transition between a low-coverage disordered phase and a c(2 x 2) phase with 
= 1/2 occurs at \i c fa 180 meV. The order parameter for the c(2 x 2) phase is |©s|. These 
isotherms are obtained from 10,000 independent samples for each value of p,. The insets show typical 
equilibrium configurations in the disordered phase at ft, — +100 meV (left) and in the ordered phase 
at +400 meV (right). 
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Figure 3: (a) Simulated CVs for L = 128 and 256 at various potential-sweep rates p. The curves 
for nonzero p give p~ 1 d<3/dt, while for p = we show dO/dp, obtained by differentiating the sim- 
ulated equilibrium coverage isotherm in Fig. |[ (b) Peak separations for various sweep rates. For 
experimental comparison, the units are given as mV and are obtained with 7 = —0.73. 
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Figure 4: Coverage, 6, and order parameter, |6s|, shown vs time for sudden disorder-to-order 
potential step at room temperature with L=256. Step from fix = —200 meV to p,2 — +600 meV, 
averaged over 10 independent runs. The inset shows the correlation length for the c(2 x 2) order 
parameter vs t 1 ^ 2 . 




Figure 5: Coverage, 0, and order parameter, |0s|, shown vs time for sudden order-to-disorder 
potential step at room temperature with L=256. The four dynamic regimes discussed in the text 
are labeled and separated by vertical bars. Step from p,\ — +600 meV to ^2 = —200 meV, averaged 
over 3 independent runs. 
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